Disorder-induced tail states in a gapped bilayer graphene 
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The instanton approach to the in-gap fluctuation states is applied to the spectrum of biased bilayer 
graphene. It is shown that the density of states falls off with energy measured from the band-edge 
as v{t) oc exp(— \t/tt\ ), where the characteristic tail energy, et, scales with the concentration of 
impurities, n„ asn/ . While the bare energy spectrum is characterized by two energies: the bias- 

1 /% 2/3 

induced gap, V, and interlayer tunneling, tj_, the tail, e t , contains a single combination V ' i ± . 
We show that the above expression for v(e) in the tail actually applies all the way down to the 
mid-gap. 

PACS numbers: 71.55.Jv, 71.23.-k, 73.20.Hb, 73.21.Ac 
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I. INTRODUCTION 

Several experimental studies of electronic proper- 
ties of graphene bilayers were recently reported in the 
literaturei2^£&l£. While experiments Refs. [Ull! 
were carried out on unbiased, and thus gapless 9 bilayers, 
the focus of the papers Refs. @@H was the fact that a tun- 
able gap emerges in the energy spectrum of bilayer upon 
applying an interlayer bias&i^. Various consequences 
of the o pening o f the gap were studied theoretically in 
Refs. [Tl|lillljllilil[iltilillilllll One of these con- 



sequences is that biased bilayer responds to disorder as a 
"normal" semiconductor, i.e., impurities give rise to the 
tails of the density of states which extend into the gap 
from the bottom of conduction and from the top of the 
valence band. Such in-gap localized states are especially 
relevant to the experiment Ref. [8|, where the inelastic 
transport over these states has been observed. This raises 
a theoretical question about the shape of in-gap fluctua- 
tion tails in bilayers and their dependence on the disorder 
strength. The only paper on biased graphene bilayers 
with impurities that we are aware of is Ref. Il3l . This 
paper studies not the tails, but rather disorder-induced 
smearing of the band-edges. Also, the numerical results 
for the density of states in the gap region are presented 
in Ref. |_13| for particular values of impurity concentration, 
rii, so that the general dependence of the magnitude of 
smearing on rii, as well as on the interlayer bias, V, was 
not established. 

In the present paper we study analytically the den- 
sity of disorder-induced localized in-gap states in bilayer 
graphene. The reason why classical result s 23 i 24 ' 25 for 
the fluctuation tails do not directly apply to this situ- 
ation, is a peculiar structure of the bare energy spec- 
trum shown in Fig. 1. In particular, the minimum (max- 
imum) of the electron (hole) dispersion is located at fi- 
nite momentum 10 , pq. Also, at energies of the order of 
the gap, the dispersion law is not quadratic, but rather 
e(p) oc p A . These two features naturally define two 
regimes of disorder-induced broadening of the density of 
states: 

(i) Weak disorder. The magnitude of smearing in this 
regime is smaller than the depth of the minimum in 




FIG. 1: Energy spectrum of a biased bilayer graphene has a 
loop of minima of depth e m Eq. (|10|l at |p| = po, Eq. (|6]). At 
small bias, V < t±, the gap is smaller than the distance to 
the next subbands. 



Fig. 1. As a result, the states responsible for the smear- 
ing have momenta close to pg. This fact, as was pointed 
out in Ref. [13|, facilitates localization of electrons (holes) 
by weak impur ities. Earlier this observation was made in 
Refs. H3|28|l29|. 

(u) Strong disorder. The mexican-hat-structure in 
Fig. 1 is completely smeared. In this regime, the in-gap 
states are formed due to trapping of electrons (holes) with 
quartic dispersion by certain disorder configurations. 

Below we demonstrate than in both regimes the mag- 
nitude of smearing, e t , is proportional to the combination 
n i V 1 ^- Clearly, for the gap to be resolved, the applied 
bias must exceed et- This suggests that the threshold 
bias, for which V > et is proportional to the first power 
of ni. We find the shape of the density of states near the 
band-edges for two regimes by extending the instanton 
approach of Refs. [23ll24l to the spectrum of a biased bi- 
layer graphene. We restrict consideration to the case of 
short-range impurity potential, w(r), so that the corre- 
lator of the disorder potential is 



(U(r)U(r'))= 7 8(r-r'), 



(1) 



with 



where 



7 = n 'i 



drw(r) 



(2) 



We also assume that the bias is smaller than the inter- 
layer tunneling constant**^!, t±. 



16tf 



(10) 



is the depth of the minimum, Fig. 1, and the dimension- 
less function v(z) is defined as 



II. BARE DENSITY OF STATES 

Due to interlayer hopping, the spectrum of the bilayer 
graphene becomes parabolic 9 -*^, e(p) = ±c 2 p 2 /t±, where 
c is the Dirac velocity in graphene. For "small" momenta, 
cp <C t±, the gap opens upon applying the bias, V, be- 
tween the layers. For V < t±, the low-energy Hamilto- 
nian of the bilayer graphene can be reduce d 11 ' 22 to the 
2x2 matrix 



H = 



(H 1 



v 



t± 



C 2 (Px+JPy) 2 \ 



2 2 
1 _ C P 



(3) 



which yields the spectrum 



e 2 (p) 



V 2 



9 2 \ 2 

c-p 1 



4 4 
C p 



(4) 



It is seen that in addition to opening the gap, finite bias, 
V, leads to negative effective mass at small momenta. 
Disorder affects the energy domains \e(p) ± V/2\ <C V, 
close to the band edges. In these domains the spectrum 
can be further simplified to 



e ± (p) = ± 



where 



V C 4 ( 2 



V 
P »=2c 



Po 



(5) 



(6) 



is the minimum position, Fig. 1. In Eq. ([5]) we also took 
into account that V < (i. The cigenfunctions, corre- 
sponding to the two branches, have the form 



X+W 



„2i4> p 



vu 



XpO) 



c *Pr 



Vlj 



(J) 



where P = tan (py/p x ) is the azimuthal angle of the 
vector p. The minimal energy for the branch e + (p) is 

V V 3 
A — -164- < 8 > 

Shifting the energy scale origin to A, for the bare density 
of states^ 5 -, vq (e), we will have 



^o( e ) = \y#) p ( e / e ™)' 



(9) 




Single-scale behavior of v(z) ensures that the magnitude 
of disorder-induced smearing of the band edges is defined 
by a single parameter for both weak- and strong-disorder 
regimes. 



III. WEAK DISORDER 

In the vicinity of the minimum at p = p$, the disper- 
sion law simplifie d 13 ' 15 ' 19 to 



6 ± (p)=± 



(IPl-Po) 2 
2m 



ti 



2Vc 2 ' 



(12) 



This expansion applies in the domain, (|p| — po) <C V/c, 
where e(p) — e(po) is smaller than the energy distance 
V 3 /16£ 2 between the minimum and maximum in Fig. 1. 
The eigenfunctions Eq. ([7]) simplify to 



xt(r) 



= 2i0 p v 



V 
'it. 



X P (r) 



«pr 



V 
it. 



(13) 



One-dimensional character of the spectrum is reflected in 
the e -1 / 2 behavior of the density of states Eq. ([9]) near 
the band-edge. The magnitude, et, of disorder- induced 
broadening can be estimated from the following reason- 
ing. The inverse scattering rate due to the disorder, cal- 
culated from the golden rule, is given by 



_l_ 7 ( ,_ -yt±v 1/2 

r(e) tt 0{> 2 7 rc 2 e 1 / 2 ' 



(14) 



Then the estimate for et emerges upon equating l/r(e t ) 
to e t , yielding 



e* 






2/3 



(15) 



Note that the above consideration equally applies for 
weak disorder, et < £m, and strong disorder, et > e m . 
As follows from Eq. (|15|) . the weak-disorder regime real- 
izes for 



7< c z 



(16) 



The fact that the fluctuation states in weak-disorder 
regime are composed essentially from the free states with 




FIG. 2: (Color online) Dimensionless density of states near 
the band edge is plotted from Eq. (|22|l versus dimensionless 
energy x = e/e t . Green line is the high-energy asymptote, 

vo(x) = Ci/x) 1 ' 2 . 



momenta, p, close to po, allows to obtain an asymp- 
totically exact solution for the density of states in this 
regime. This was demonstrated in Ref. [28|, where the 
spectrum Eq. (j 1 2[) emerged as a result of the spin-orbit 
coupling. The difference between the wave functions 
Eq. (|13p and spin-orbit wave functions affects only nu- 
merical factors in the final result. Still, for completeness, 
we will briefly sketch the calculation from Ref. [28|, using 
our notations. 

A drastic simplification coming from the condition 
et <C £ m is that the interference between two scattering 
amplitudes p — > p' and p — > pi — » P2 — > p' is sup- 
pressed. This is because the momenta p, pi, P2, which 
are close to po i n absolute value, are restricted in their 
mutual directions by the condition (|pi + P2 — p| —po) ~ 
Poi^t/cm) 1 ^ 2 - This condition limits the angles between 
the momenta to (et/e m ) • As a result, the non-random- 
phase-approximation (non-RPA) diagrams in the self- 
energy, Sp(e), are parametrically, in (fit/em) 1 > smaller 
than corresponding RPA diagrams. In other words, 
the RPA becomes asymptotically exact in weak-disorder 
regime. 

Within the RPA, the density of states is given by 
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1 
= —Irr 


p 


Ixp +, I 2 

-e(p)-S p (e)' 


where the electron self- 


energy 


satisfies 


f m v (rr.\ — r-j 


f 


dpi 


\(x P M x P V)\ 2 



(17) 



) 2 S-e( Pl )-E Pl (i?)- 



(18) 



The fact that IniE p does not depend on p allows to 
express the solution of Eq. (fT8|) in the form 



ImS(e) 



&* \ e t ) 



(19) 



where the 


energy e 


is 


defined as 




e = E 




(t- 


V 3 

164 



ReTi\ , 



(20) 



and the dimensionless function A(x) satisfies the alge- 
braic equation 



A(x) = 



VM*) 2 



A(x) 2 + x 2 
The solution of this equation has the form 



A(x) 



3 i/2 2 i/3 (33/2 + ^7 + 4x 3 ) 



1/3 



(3 3 / 2 + V27 + 4x 3 ) 



2/3 



2 2 / 3 ir 



(21) 



1/2 



(22) 



The density of states per spin and per valley can be ex- 
pressed via A(x) as 



1/(6) 



1 

7T7 



Im'E(e) 



1 



VV^j. /2 5 / 3 e 



2 5 / 3 7T 1/2 o 



e< 



(23) 



It is plotted in Fig. 2 together with the bare density of 
states. 

Sharp low boundary of u(e) in Fig. 2 is an artifact 
of the RPA. In reality it is smeared within the energy 
interval 



et 



3/2 

1/2 



474 
Vc 2 



<e t , 



(24) 



which is much smaller than et . This smearing comes from 
non-RPA diagrams. Concerning the deep tail of v(e), it 
is close to a simple exponent, namely 



v(t) (x exp 



4|e| 



e t ln(e m /|e|) 



(25) 



The reason why this tail can be found analytically is 
again the fact that the wave functions of the fluctua- 
tion states have two spatial scales; they oscillate rapidly 
with period p^ and decay at much larger distances as 
exp(— ^/2m|e| r). The tail states are very similar to those 
found in Ref. [2f|. The key steps of derivation of Eq. ([25]) 
are outlined in the Appendix. 



IV. STRONG DISORDER 

We now turn to the case of a strong disorder when 
et > e m . The shape of the tail density of states in this 
case can be established from the following qualitative 
consideration. The probability density to find a fluctua- 
tion U(r) is given by 



P{C/(r)} = exp 



-i/'/w-m 



(26) 



In order to create a localized level with binding energy, 
e, the magnitude of the fluctuation must exceed |e|, while 
the size cannot be smaller than the de Broglie wave 
length, r e , of a free electron with energy, e, i.e., 



1 

r, = — = 



Vt 2 



1/4 



(27) 



where the last identity follows from the dispersion law 
e(p) = cp/Vt 2 ,. Now the integral J drU 2 (r) can be 
estimated as e 2 r 2 . Substituting this estimate into Eq. 
(|2"rJj) and using Eq. fl27|), we get 



v(e) ex exp — 



3/2 



(28) 



The remaining task is to establish the numerical coef- 
ficient in the exponent Eq. (j2"8|) with the help of the 
instanton approac h 23 ' 24 i 25 . Within this approach, one 
should solve the Schrodinger equation with potential, 
U(r), which yields an eigenvalue, E (here we measure 
energy from the gap center). Then U(r) is determined 
from the condition that J dr U 2 (r) in the exponent of 
Eq. (f2"o| is minimal. This restriction is convention- 
ally incorporated by adding to JdrU 2 (r) the energy, 
E[%] = (ifr\(H+U)\&) with Lagrange multiplier, A. Then 
minimization of 



\($\(H + U)\V)- j drU 2 (r) 



(29) 



with respect to U yields U = — ^l^l 2 . 

At this point the following remark is in order. Con- 
ventionally, upon substituting the found U(r) back into 
the Schrodinger equation, the sign of A is chosen from 
the condition that potential U(r) is attractive. However, 
in our case of two symmetric bands, Fig. 1, the potential 
which is attractive for electrons, is repulsive for holes, and 
vice versa. It turns out that choosing A = 2 corresponds 
to v(E) which falls off from E = V/2 all the way down 
to E = — V/2, Correspondingly, choosing A = — 2 leads 
to the tail v(—E) which grows from E — —V/2 towards 
the bottom of conduction band E = V/2. Therefore, at 
E = 0, we have two different solutions, with U(r) and 
— U(r), for which the "electron" and "hole" components 
of eigenfunction \& are related as (ip e ,tph) *-* (— V'LV'e)- 
Thus, in view of exponential character of v(E), it is suf- 
ficient to set A = 2 and consider only positive energies, 
E > 0. The nonlinear instanton equation reads 



^Ve " Y (d x + idyf i> h - tP e (l^l 2 + \tP h \ 2 ) 



Eip e ; 

E^ h . 
(30) 



Then, with exponential accuracy, 
V{U(t)}, where 



WO 



have v(E) 



V{U(r)} = cxp 



-i/./r^lrj 



(31) 



t 
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FIG. 3: (Color online) The function </>(s) is plotted from 
Eq. (35) . 



exp 



~Jdr(\A 



|Vv* 



2\ 2 



is the probability of realization of U(r). 

Consider first the energies close to the bottom of con- 
duction band, — e = (¥ — E) <C V. In this limit, the 
second equation in the system Eq. (|30|) can be simplified 
as 



iph = 



(d x - id y ) ip e 



Vt 



(32) 



Substituting Eq. ([32)1 into the first equation Eq. (|30|) . 
and performing rescaling 



r = c{Vtl\e\)- l ' i p, ^ e = le/AI 1 / 2 /^), 



(33) 



we arrive to the following dimensionless instanton equa- 
tion 

A 2 p f( P ) + f(p)-f(p) 3 =0, (34) 

while the expression for v(e) takes the form 

3/2" 



v{e) ex exp 



h 

2 



(35) 



where et is defined by Eq. (fT5jl and I4 = J dpf 4 (p). 

Eq. (|30[) has solutions for arbitrary angular momen- 
tum, M. However, leading contribution to the density of 
states comes from azimuthally symmetric solution, f(p). 
Then the hole component Eq. (|32p of the wave function 
corresponds to the momentum M. — 2, namely 



^.-(«:-i«.)/w. (36) 



iph 



A peculiar feature of the non-linear equation Eq. 
is that it contains A 2 instead of a usual Laplace opera- 
tor, A p . This is a direct consequence of the dispersion 



CONCLUDING REMARKS 




FIG. 4: (Color online) Electron and hole components of the 
wave function, obtained from variational solution of Eq. (|30[) 
for E = 0, are plotted versus dimensionless distance, p, de- 
fined by Eq. 



law e(p) a p 4 . As a result, the average "kinetic" en- 
ergy, J 2 = f dp[A p f(p)] 2 , "potential" energy I4 and the 
integral I2 = J dp[f(p)] 2 are related as 1:4:3, unlike 
the relation 1:2:3 for conventional polaroid. We solve 
Eq. (|34p by employing the variational approach on the 
class of trial functions f(p) — Cexp[— (p/p ) s ]. Mini- 
mization of the corresponding functional 



${/} 



dp 



(A P /) 5 



1 



f 



(37) 



with respect to C and po can be easily performed analyt- 
ically. The resulting s-dependence of $ {/} has the form 



0(a) 



_s(s-l) + 2 (1-2/a) 
23-4/s sin7r(l-2/s) 



1/2 



(38) 



This combination has a well-pronounced minimum atsw 
1.963 see Fig. 3. On the other hand, the "hole" wave 
function Eq. (|36|) is non-singular at p — > if s > 2. Thus 
it is reasonable to adopt s — 2, yielding C = 1.63 and 
po = 2.21. This allows to specify the numerical factor 
in the exponent of density of states in Eq. (|35|) . namely 
h = 27.36. 

The result Eq. (|35[) applies for e t < e < I'. A relevant 
question for inelastic transport is the density of states at 
the gap center. For qualitative estimate it is sufficient 
to substitute e = — V/2 into the Eq. (J35J), which recovers 
Eq. (??). To establish the numerical coefficient in the 
exponent more accurately, we found variational solution 
of the system Eq. (|30|) for E = 0. It turns out that the 
coefficients C and po in f(p) assume the values 1.6 and 
1.96 respectively, and the numerical coefficient in Eq. ([35]) 
is I 4 = 23.22. We conclude that the expression Eq. ([33]) 
for the tail of the density of states is essentially valid 
down to the gap center. 



1. Prefactor. Eq. (|35|) describes the tail of the density of 
states with exponential accuracy. To restore the dimen- 
sionality, it is natural to multiply Eq. (|3l))) by vo(e t ), since 
the smearing of the band-edge is ~ Ft- However, since 
Eq. (j3"5|) describes the deep tail, the prefactor contains 
an additional power of a dimensionless ratio (|e|/et). To 
establish this power, one has to follow the procedure of 
calculating the prefactor in the functional integra l 30 ' 31 . 
Within this procedure, the origin of a prefactor is the 
fact that the center of the instanton can be shifted in 
the plane along both axes. These shifts correspond to 
the so-called zero modes. Each zero mode contributes a 
factor (lel/et) 1 / 4 . The power 1/4 reflects the size of the 
instanton fluctuation, r e oc lei" 1 / 4 , see Eq. (|2T|) . Overall, 
within a numerical coefficient, the final expression for the 
density of states in the tail reads 



< ^ Vl/H ± > ,1/2 



c 2 t t 





F 


3/2" 


-11.6 








Q 





(39) 



Note that, unlike the case of parabolic spectrum^ ., for 
f(p) oc p 4 , the prefactor does not diverge. This is because 
the second-order shift of the band-edge oc 7 J cfei/ ( e )A 
converges at large e for i'o(e) oc e" 1 / 2 . 
2. Relation to the scattering time. Short-range disorder 
is characterized by a single parameter, 7. For compar- 
ison with experiment, this parameter can be related to 
electron scattering time, t>, in the case when the gate 
voltage places the Fermi level, Ep, well above the smear- 
ing, Ft, of the band-edge. Expressing 7 from Eq. (fl4)l . 
and substituting into Eq. (|15[) . we obtain 



et 



/ 4tt 2 ^ 



1/3 



47T 4 n 2 c 



1/3 



(40) 



In the second identity we had expressed E F via electron 
density, n e . 

3. The role of intervalley scattering. Our assump- 
tion that disorder is short-ranged requires that the ra- 
dius of the impurity potential, w(r), is smaller than 
the wavelength of electron with energy ~ V, which is 
~ c/ (Vt±) . This length is much larger than the lat- 
tice constant, and we neglected the intervalley scattering. 
If the radius of w(r) is comparable to the lattice constant, 
the intervalley scattering becomes as efficient as intraval- 
ley scattering. This would not only lift the valley de- 
generacy of the fluctuation states but also result in their 
azimuthal asymmetry, much like in the case of degener- 
ate valence band considered in Ref. [33|. The consequence 
of this asymmetry is the change of the numerical factor 
in the exponent of Eq. (133)) . 

4. Dependence on impurity concentration. Our main 
result Eq. (|39|) applies in the limit of strong disorder when 
7 > c 2 (V/t±) 4: . On the other hand, we assumed that the 
gap is not washed out completely by the disorder. Then 
the upper limit on 7 can be found by setting e = V/2 



in Eq. (|39j) equating the exponent to 1. Finally, it is 
convenient to present the domain of validity of Eq. 



7 V 



(41) 



It follows from the second condition that the minimal 
V = V c , at which the gap effectively opens, is propor- 
tional to the impurity concentration, as it was stated in 
the Introduction. We can also rewrite the above con- 
dition in terms of dimcnsionless conductance when the 
Fermi level is in the conduction band 



E f t f > 1.5 



El 

V 



3/2 



(42) 



Finally, we address the case of low impurity concentration 
when the first condition Eq. (|4"Tj) is violated. 

An isolated impurity with potential, w(r), creates a 
localized state with binding energ y 27 i 28 ' 29 



impurity concentration, i.e., that the levels Eq. (|43)l are 
not formed at n* satisfying the condition Eq. (j45|) . In 
conclusion, we rewrite for completeness the condition of 
strong disorder in terms of impurity concentration and 
binding energy of an individual impurity 






< 



TT 2 p 2 



<4.1- 



t± 



(47) 



We see that the smaller is the gap, the broader is the 
interval Eq. ([47]) . 
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APPENDIX A 



e 6 = 2n mp / drw{r) J (p r) 



(43) 



For the short-range potential, w(r), the Bessel function 
in the integrand can be set to 1. Substituting p a and m 
from Eqs. ([6]) and (fl"2"l) , we obtain 



q< 



TT 2 Vt 2 A 

Ac 2 



drw(r) 



(44) 



We see that, for bilayer graphene, the binding energy 
is proportional to the gap 13 . The wave function of this 
localized state not only falls off exponentially with dis- 
tance, r, from the impurity, but also oscillates as Jo(por). 
It is clear that when the average distance between the im- 



purities 



-1/2 



becomes smaller than p the impurity 



band merges with the conduction (valence) band, 
markably, the criterion 

m > pi 



Re- 



(45) 



also follows from a very different reasoning. The expres- 
sion Eq. (|2~4"|) for the tail energy, e t , in the weak-disorder 
regime can be rewritten as 



e* 



An t t\ 
Vc 2 



drw(r) 



(46) 



Comparing Eq. (|4"6"|) to Eq. (|4"3"| . we find that the ra- 
tio et/eb is ~ nt/Po- Thus the condition Eq. (|45|) that 
neighboring localized states overlap ensures that the im- 
purity band transforms into the tail of the conduction 
(valence) band. Note also, that in treating the weak- 
disorder regime we assumed that the correlator of the 
disorder potential is given by Eq. {T]). By making this 
assumption we already implied that disorder is not due 
to individual impurities but rather due to fluctuations of 



Similar to Eq. pTjl . i/(e) in the tail is given by 

v(e)<xexp(-±-jdr\cp(v)\ 4 y (Al) 

where the two-component function ip(r) satisfies the 
equation 



H<p(r)-\<p(r)\ 2 <p(r) = e<p(r), 



(A2) 



with H being the free Hamiltonian with the spectrum 
Eq. (fT2|) and the eigenfunctions Eq. (|13[) . Searching the 
solution in the form 

^(r) = ydpB(p) X +(r), (A3) 

we arrive to the following integral equation for B(p): 



B(p)[ e +(p)- e ] = 



w/*(/JL* ,DW ) 



x IXp XpiJIXpaXps 



(A4) 



Assuming that B(p) depends only on the absolute value 
of p, we easily perform the angular integration in (|A4[) . 
Using explicit forms of the wave function scalar products, 
we obtain 

B(p)[e + (p) - e] = n 2 fdrr I f JJ d Pl B(p z ) 

\ i=l,2,3 I 



x J2(pr)J 2 (pir)J 2 (p2r)J2(P3r), 



(A5) 



where J2(x) is the Bessel function of the second order. 
The product of J2{Pif) manifests the difference of Eq. 
[5| from the corresponding equation in Ref. |28| . 



The principal step in solving Eq. (|A5|) is setting all 
momenta in the right hand side equal to Pq. Then the 
integral over r yields ^ln(e m /|e|) so that Eq. (|A5j) re- 
duces to 



B(P) 



(p-Pa) 2 
2m 



= 2p ln[^ 



dp'B{p') 



(A6) 



This equation has an obvious solution of the form 

(P> ( P -po) 2 /2m+\e\- 
Substituting Eq. (|A"7|) into Eqs. (|A"6J) and (|A"3l) . we find 



(A7) 



for constant j3 the value 



-3/4 



P = ^km^' 2 (^y ^'WW), (A8) 

and for ip(r) the form 



tp(r) = 27r 2 /3p 



1/2 / J 2 (por) 







(A9) 



which is valid for r < (2m|e|) l / 2 . Finally, Eq. 
emerges upon substituting Eq. (|A9|) into Eq. (|A1|) . 
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